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We present a framework for phenomenological lattice QCD calculations which makes use of a 
tree level Symanzink improved action for gluons and stout-link Wilson fermions. We give details of 
our efficient HMC/RIfMC algorithm and present a scaling study of the low- lying Nf — 3 baryon 
spectrum. We find a scaling region that extends to a < 0.16 fm and conclude that our action and 
' algorithm are suitable for large scale phenomenological investigations of Nf = 2 + 1 QCD. We expect 

this conclusion to hold for other comparable actions. 
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I. INTRODUCTION 



Over the last decade, it has become clear that smeared-link fcrmion actions^ offer substantial technical advantages 
over their thin-link counterparts. The idea of damping unphysical UV fluctuations by replacing elementary links with 
[ a weighted sum of paths was first introduced in the framework of pure gauge theory T|. It was later recognized 
^ ■ that the chiral properties of clover fermions 0] can be substantially improved by replacing the thin links in the 
Oh' covariant derivative of the fermion operator with their smeared counterparts [11 . From a Symanzik point of view, this 
^ . replacement amounts to adding ultralocal irrelevant terms to the fermion action, as long as the smearing prescription 
(parameter, iteration number) stays fixed as a function of bare coupling. In this way it is guaranteed that the 
continuum limit is unchanged. 

CN ■ In the context of quenched QCD, the advantages of smeared clover fermions are well established [1, S B 0, Hi ■ The 
K*" I theoretically leading 0{asa) contributions are, in practice, absent and the extrapolation to the continuum appears to 
I be dominated by 0{a^) cut-off effects. In particular, the tamed UV fluctuations result in improved chiral symmetry 
\ properties. Furthermore, the smearing significantly reduces the contributions of unphysical tadpoles; renormalization 
. constants are generally closer to their tree level values, and csw is not far from 1 at typical lattice spacings. 
• ■ Given this experience, it is reasonable to expect that also dynamical clover fermions will benefit from link smearing. 
There, the non-differentiable nature of the back projection step of the smeared link onto the gauge group, which is 
' usually performed, for instance, when using APE smearing [l|, may pose problems for the molecular dynamics update. 
[ An early suggestion was to use "stout" links 0] to define fermions which can be simulated with the Hybrid Monte 
' Carlo (HMC) algorithm [l^. Further particulars of the HMC force with UV-filtered actions have been worked out 
. !^ ' in Recently, several alternative smearing methods suitable for dynamical simulations have been proposed [l^ . 

^ ' The efficiency of link smearing results from the fact that it leaves the structure of the fcrmionic operator entirely 
unchanged. Smeared clover fermions still have exclusively nearest-neighbour couplings. The damping of unphysical UV 
modes is achieved exclusively by a modified - but still ultralocal - coupling to the gluonic background. This modification 
of the fermionic action is continuum irrelevant, but at a given finite cutoff one generally expects observables with weaker 
coupling to unphysical UV modes to be closer to their continuum limit values, resulting in overall improved scaling. 
In the present paper we investigate this issue by performing a scaling study with Nf = 3, stout-link clover fermions. 
Although other smearing methods are presently known, we opt for the standard stout-link prescription because it is 
widely used and will share features, such as an enlarged scaling region, with other comparable prescriptions. 

The size of its scaling region is one of the most important criteria to assess the suitability of a given action for 
phenomenological purposes. The onset of scaling, together with the power of the lattice spacing against which results 
need to be plotted to show a linear dependence, determines the finest lattice spacing needed to reliably extrapolate 



* CPT is "UMR 6207 du CNRS et des universites d'Aix-Marseille I, d'Aix-Marsoillc II ot du Sud Toulon- Var, affilioo a la FRUMAM". 
^ In the literature they are also referred to as "UV-filtered" or "fat-link" actions. Likewise, actions in which the covariant derivative 
involves the original gauge links arc sometimes called "thin-link" actions. 
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FIG. 1: Performance of CG in double precision (squares) compared to a mixed precision variant of CG (circles). Data are 
from an iV/ = 2 + 1 run on a 32^ x 64 lattice at /3 = 3.57 with amJ^J^° - 0.0077 and ani^'~^^'^ ~ 0.049 corresponding to 
~ 250 MeV. 



to the continuum, and hence the overall cost in terms of CPU time. Our main result is that smeared clover fermions 
do indeed show very nice scaling properties up to at least 0.16 fm lattice spacing. Moreover, the link-smearing seems 
to eliminate known pathologies that unfiltered actions may show in a dynamical setting [13, [ill . 

The results presented in this paper are obtained using a tree-level Symanzik improved gauge action and 
six-step, stout-smeared clover fermions with a clover coefficient taken at its tree-level value cgw = 1 (though a 
perturbative 0, [Ml or non-perturbative [H, determination is feasible). Note that also in the clover term 
is built from the same set of stout links. This choice allows for efficient simulation while delivering good scaling 
properties, as demonstrated below. Moreover, dedicated studies in quenched QCD have shown that the dependence 
of observables on smearing is quite mild (see e.g. [I,[l3|) and the exploratory studies of e.g. [HI, [13, H, [3 suggest 
that this behavior persists in the the full theory. Thus, our choice involves no fine-tuning and we expect our results 
to hold for actions which involve comparable amounts of smearing. 

In our scaling study we choose Nf — 3 for simplicity, creating an artificial world with degenerate u, d and s quarks. 
We will denote the pseudoscalar and vector mesons by tt and p respectively. Our goal is to perform continuum extrap- 
olations along three distinct lines of constant "physical" quark masses, characterized by M^/M^ = 0.60, 0.64 and 0.68. 
Since we do not aim in this paper at phenomenologically relevant computations and instead would like to test the 
extent of the scaling regime, we deliberately choose these rather large masses. With Ma for standard hadrons close 
to one, cut-off effects with inferior actions will be large. In all our runs M^L is kept fixed, at values larger than four, 
to avoid finite volume effects. Our goal is to simulate at several values of the gauge coupling and fixed M^/Mp and 
A'/^L, and to determine the scaling of the baryon octet and decuplet masses, Mn and A'/a- 

The remainder of the article is organized as follows. In Sec. 2 details of the action and our algorithm are given. 
Sees. 3 and 4 are devoted to tests which provide clear evidence for the absence of bulk phase transitions in our 
simulations. In Sec. 5 we show that our action is ergodic with respect to topology. Sec. 6 then contains a detailed 
scaling study of the nucleon and delta masses. We conclude with a short summary and outlook. 

II. ACTION AND ALGORITHMS 
A. Action 

The explicit form of our gauge and fermion action in terms of the thin (Un.fj.) and smeared (Vn^p,) gauge links is as 
follows: 

C _ cSym , qSW 
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^Sym ^ ^ ^ Re Tr (1 - C/piaq) + |- ^ Re Tr (1 - (7,ect)] (1) 

plaq icct 
n 

with the standard Wilson action 5*1^. The parameters csw: cq and ci set to their tree level values: 

csw = 1, ci = -1/12, Co = 1 - 8ci = 5/3. 

Both the hopping part and the clover improvement term in the fermion action use six-step stout-smeared links @ 
Vn,ii = KiS- Those are constructed from the thin links Un,^L = according to 

5'(") — i(p(")^(")t _ -^/(")p(n)t^ _ iReTr(r''"''p^^"''''' — y(")r^"^''') (2) 
2 6 

n,fjL / J n,v n-t-i/,/1 n-t-/x,i^ 



The stout smearing parameter is chosen to be p = 0.11, which is a rather conservative choice [1, Q corresponding to 
an aAPE = 0.48 with respect to the average plaquette [l2|. In S^'^ only the unfiltered links are used. As detailed in 
the Introduction, this action is ultralocal in both the quark and gauge sector. 



B. Simulation algorithm 



We start with the description of our Nf = 2-1-1 algorithm. Two flavors are implemented via the Hybrid Monte 
Carlo (HMC) algorithm [l3|, the third using the Rational Hybrid Monte Carlo (RHMC) algorithm [13, Hi]. We 
employ even/odd preconditioning (25| to speed up the fermion matrix inversions. The generic HMC algorithm suffers 
from critical slowing down in the light- quark regime. To treat this problem, we combine several improvements over 
the generic algorithm (see also [26l |27|): 

• Multiple time-scale integration: not all force contributions in the molecular dynamics (MD) part of the HMC 
algorithm require the same amount of computational resources. Using multiple time-scale integration ("Sexton- 
Weingarten integration scheme" ) [28^ , it is possible to put each part of the MD on a different time scale according 
to its relative contribution to the total force, thus reducing the computational costs of the MD. 

• Mass preconditioning: the pseudofermion force is used within the MD to include the effects of dynamical 
fermions. Through mass preconditioning, the UV part of the force can be split off and treated separately [29j . 
which helps reducing the fluctuations in the force. The second important benefit of mass preconditioning 
appears when combined with the multiple timescale integration scheme [26l . [27j : the more expensive infrared 
part contributes less to the total force and can be integrated with larger time steps. 

• RHMC: the third, unpaired quark flavor is implemented through the RHMC [H, [IJ algorithm. This algorithm 
makes use of the fact that the single fermion action can be written as {M^ M)~^/'^^, where the inverse square 
root can in turn be approximated by a rational approximation and be efficiently calculated with a multi-shift 
solver. The RHMC is highly efficient in simulating a single quark flavor. It can also be combined with the 
multiple timescale integration scheme. 

• Omelyan integrator: the MD integration within the generic HMC algorithms uses the leapfrog integration 
scheme. It proceeds by first integrating one half step in position space followed by a full step update of the 
conjugate momenta and finally another half step in position space. The Omelyan integrator adds a small 
momentum update (reduced by A w 0.193) before and after the leapfrog step and shortens the original leapfrog 
momentum update in by a factor (1 — 2A). This scheme improves the MD energy conservation by about one 
order of magnitude for a factor ~ 2 increase in computational cost. The use of a correspondingly larger step 
size then results in a net gain of about 50% [30] . 

We use this algorithm also for our Nf = 3 scaling study with mnMC = 'tirhmc- 
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C. Inversion algorithms 

The most time consmxiing part, both in the valence and the sea sector, is the (approximate) fermion matrix inversion 
by means of a Hnear solver. These calculations generally require double precision accuracy. This is due to the fact 
that, in order to maintain reversibility, the MD part of the algorithm has to be performed in double precision. Double 
precision accuracy is also required in valence calculations at small quark masses, owing to the large condition numbers 
involved. However, this does not imply that each fermion matrix multiplication needs to be done in double precision. 
In the valence sector we need to solve 

Dx^b (3) 

(with D in our case being the stout-link clover Dirac operator) to construct the correlators. To calculate the fermionic 
force in the MD part of the algorithm we need to solve 

ZjtDa; = b. (4) 

In both cases it is possible to use a single precision version of D within mixed precision solvers to accelerate the 
inversion. There is basically no penalty in terms of the iteration count: we find that the increase in the number of 
matrix multiplications is well below 10%. 

A simple and reasonably efficient way to construct a mixed precision solver is to use the standard "iterative 
refinement" technique, which amounts to repeatedly using a single precision solver. In this scheme, only the (outer) 
residuals and global sums are calculated in double precision; the inversion is performed with single precision accuracy. 
The single precision inversion typically uses the same algorithm that would be used for a full double precision inversion, 
such as BiCGstab to solve ([3]) or CG for Q. With A — _D or A = D'^D referring to the forward multiplication routine 
in double precision, a the single precision counterpart and e the desired final double precision accuracy, the complete 
procedure reads: 

1. Compute Ti = b — Axi 

2. If \ri\ < e\b\, exit 

3. Solve ati — in single precision to an accuracy e', with U denoting the solution. 

4. Update Xi+i = Xi + ti 

5. Goto 1 

With Si — Vi — Ati and S = \si\/\ri\ w e' < 1, we have 

\n+i\ = \b- Axi+i\ = \b-Ax, - AU\ = \b- Axi - r, + Si| = \si\ = S\ri\ < \ri\ . (5) 

Thus, as long as the single precision inversion does not fail, the method will converge. Since many single precision 
matrix multiplications are needed to compute ti, compared to just one double precision multiplication with A in the 
outer iteration, the whole solver is dominated by the single precision matrix multiplication performance, resulting in 
a significant speedup over a full double precision inversion (see Fig. [T]). 

III. SPECTRAL GAP 

In quenched QCD, the (unsmeared) clover fermion operator may have one or several eigenvalues close to the origin 
or with a negative real part, even for not very light quark masses. Configurations for which this is the case are referred 
to as "exceptional" . 

If one integrated the HMC trajectories exactly, any such configuration would be absent in full QCD, since an 
eigenvalue of the hermitean Wilson operator Hy/ — 75 -Dw approaching zero would induce an infinite back-driving 
force in the HMC. In practice, when the trajectories are generated with a finite step-size integrator, the near zero 
modes along a trajectory are only approximately suppressed. This may cause a breakdown of the MD evolution. It is 
therefore natural to monitor the smallest eigenvalue (in magnitude) of flw a-nd check if it is sufficiently far from the 
origin throughout the entire run. In a given ensemble this spectral gap shows a more-or-less Gaussian distribution, 
and as long as its median is several a away from zero, the simulation is deemed safe (sij . 

Since we use even-odd preconditioning, the relevant quantity to monitor is the smallest eigenvalue of the hermitean 
counterpart of the reduced operator Dj-cd = ^(^^oo ~ ^oc-^oc^-^co)' which is 75-hermitean. We include a factor 1/2 
to have its IR eigenvalues almost aligned with the low-lying eigenvalues of the full operator. For the lightest mass 
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FIG. 2: The magnitude of the smaUest eigenvalue of the preconditioned hermitean Dirac operator in units of the PCAC mass. 
At each /3 the hghtest run {AItt/Mp ~ 0.6) is shown. 
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FIG. 3: Histogram of the inverse iteration number of our Unear solver at a lighter Af^ for the lightest pseudofermion in the 
action. Results are from an TV/ = 2 + 1 run on a 48^ x 64 lattice at /3 = 3.57 with am^f -^"^ ~ 0.0056 and amf '^'^'^ ~ 0.044 
corresponding to Mtt ~ 190 MeV. 



{M-^/Mp = 0.60, of. Sect. 5) the distributions are shown in Fig. [21 with (3 ranging from 2.8 (left) to 3.76 (right). One 
can see that even for the strongest coupling, there is still a clear separation of the eigenmodes from the origin. 

For phenomcnological applications it is of course most relevant to know how this spectral gap evolves when lowering 
tlie masses of two of the three flavors. Instead of monitoring the lowest eigenvalue of 75Dredj we opted for monitoring 
the closely related quantity l/ncc: where ncG is the iteration count for the lightest pseudofermion in the action 
for our Nf — 2 + 1 runs. In Fig. [3l we plot a histogram of l/ncc for one of our lightest production runs (for 
phenomcnological studies) and find a clear gap, which provides strong evidence for the stability of the algorithm. We 
have also monitored the acceptance rate and the Hamiltonian violation AH throughout our runs and have seen no 
sign of any algorithmic problems. 



IV. SEARCH FOR POTENTIALLY METASTABLE BEHAVIOR 



In dynamical Wilson fermion simulations with small quark masses, it was reported that the system appears to 
undergo a first-order transition to an unphysical phase [3 IS^I • This was argued to mean that there is a lower bound 
on the quark mass, below which physically sensible simulations cannot be performed. Moreover, it was observed, that 
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FIG. 4: Absence of hysteresis in the average expectation value of the plaquette. Data are from an Nf =2 + 1 run on 
a 16^ X 32 lattice at /3 = 3.3 with a fixed strange quark mass am^'~'**^ — 0.0677 and the light quark mass varying between 
am^ J**^ — 0.0066 and 0.0243 in ascending (square) and descending (circles) order. The range of light quark masses corresponds 
to Mtt ~ 240 — 440 MeV. The second data set is slightly offset along the x-axis for better readability. 



1. the phenomenon occurs only with coarse lattices, 

2. gauge action improvement decreases the lower bound on the quark mass [s^, 

3. 0(a)-improved Wilson fermions together with improved gauge actions made the problem disappear for all lattice 
spacings investigated in jlS] . 

4. one level of stout smearing weakens the phenomenon [3^ . 

When discussing such phenomena, it is important to remember that a first-order phase transition can only occur 
in infinite volume. In finite volume, the metastability can be understood as an artifact of the updating algorithm: 
with an efficient algorithm, the system should eventually find the true minimum of the effective potential. Thus, for 
finite-volume simulations, the relevant question is: can the algorithm thermalize the system in a manageable number 
of updating steps? 

To investigate this issue, we have taken two 16'^ x 32 configurations, one with random links and the other, ther- 
malized in a Nf = 2-1-1 simulation at /3 = 3.3, with am^^^^ = 0.0066, corresponding to a pion mass of approx- 
imately 240 MeV, and am^^^*-^ ~ 0.0677, corresponding roughly to the physical strange quark mass. A "down- 
ward" updating sequence was then constructed from the random configuration: consecutive simulations at o.m^^'^J^'^ ~ 
0.0243,0.0173,0.0131,0.0086,0.0066, corresponding to a range of pseudoscalar masses ~ 440-240 MeV, were per- 
formed, with each simulation starting from the last configuration of the previous (larger mass) run. Similarly, an "up- 
ward" sequence of five simulations was obtained, beginning with the configuration thermalized at am^*^^^ ~ 0.0066, 
and ending with a run at am^'"/^'" — 0.0243. For each point in the two sequences, approximately 400 trajectories 
were generated, of which the first 100 were discarded when calculating the average expectation value of the plaquette. 
The resulting plaquette values, obtained during the two updating sequences, are shown in Fig. |4l No sign of hystere- 
sis is observed: the algorithm evolves the system to the correct equilibrium state in a reasonable number of steps, 
independently of the starting configuration. 

This absence of evidence for metastability, together with the good performance of our algorithm in all of our 
production runs, gives us confidence that our choice of algorithm and of action is appropriate for the range of 
parameters that we have considered so far. 
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FIG. 5: History of the unrenormalized gluonic topological charge (left) and the corresponding autocorrelation function plot 
(right), measured on our finest lattice with the smallest quark mass: (3 — 3.76, aM^ — 0.2019(20). The integrated autocor- 
relation time of gnai is approximately 2 configurations on this ensemble. A separation of one configuration corresponds to 10 
HMC/RHMC trajectories. 



V. TOPOLOGY 



In phenomenological applications, the combined choice of an action and an algorithm must allow for an adequate 
sampling of sectors of different topological charge. Quite generically, this sampling becomes more difficult as the 
continuum limit is approached. Thus, as the lattice spacing is reduced, the autocorrolation time of topological charge 
increases. This is also the case in our simulations. However, within the range of lattice spacings which we consider, 
we observe no dramatic slowing down of tunneling events. 

To determine the topological charge of our configurations, we use the naive gluonic charge definition 

9-' = E , (6) 

X 

where is the gluonic field strength tensor and the sum extends over all lattice sites. We calculate Ffj_i, at each lattice 
site as follows. After applying our smearing prescription ([2]) to the links, we average the four plaquettes emanating 
from this site and which lie in the fi-v plane. The field strength tensor is then defined as the anti-hermitian part 
of this average. The charge defined in Eq. ([6]) leads to non-integer values and must be renormalized for quantative 
studies of topology. However, such a renormalization is not necessary here since we are only interested in verifying 
the topological ergodicity of our simulations. 

The simulation-time evolution and autocorrelation of this unrenormalized topological charge are shown in Fig. [5] 
for our finest lattice and its smallest quark mass, aM^^ = 0.2019(20). The integrated autocorrelation time is around 
2 configurations. The autocorrelation decays very rapidly and is compatible with zero within the error bars after 
around 5 configurations. We can easily conclude from these two plots that there is no long-range correlation. 



VI. SCALING STUDY 



For our scaling study, we use lattices with approximately constant physical volume at five different lattice spacings. 
We opted for an Nf — 3 instead of an Nf = 2 setting in order to test the full RHMC algorithm that is also being used 
for phenomenological applications. We choose a T = 2L geometry with lattice sizes varying from L/a = 8 to L/a = 24 
and bare gauge couplings between /? = 2.8 and (3 = 3.76. We measure fermionic observables every twenty trajectories 
for L/a = 8, 10, 12 and every ten for L/a = 16, 24. For the error analysis, we use the "moving-block-bootstrap" [s^ 
technique with a binlength of two times the integrated autocorrelation time of the quantity which is measured. This 
binlength is typically around 2 for the coarsest lattices and around 8 for the finest lattices. The number of bootstrap 
samples is chosen to be 2000, because the calculated bootstrap errors saturate at ~ 1500 samples. 

At each lattice spacing we simulate a number of masses (from seven at L/a = 8 to three at L/a = 24) such that 
Mjr/Mp is between 0.60 and 0.68. As already mentioned in the Introduction, it is preferable to use these rather large 
masses for a scaling study in order to enhance possible discretization effects of order Ma. After fixing to Coulomb 
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FIG. 6: Effective masses of the pion, rho, nucleon and delta on our ensemble with L/a = 16, /3 = 3.59, ampcAC = 0.04608(12). 
The points are obtained by solving, for aMcft, the equation C{t - l)/C{t + 1) = f{aM^s{t - 1 - T/2)) / f{aMcs{t + 1 - T/2)) 
at each t, where f{x) — cosh(a;) (for tt and p) or f{x) = sinh(a;) (for and A). The horizontal lines are the masses with error 
bars obtained from correlated cosh or sinh fits to the corresponding two-point functions in the time intervals indicated by the 
length of the lines. 

gauge, we measure propagators with multiple Gaussian sources on different time slices. The source size is set to Z//4 
and is thus roughly constant in physical units. Using a Gaussian sink of the same size, the effective masses usually 
reach a plateau very quickly and we can determine a useful fitting window from it. Note that - because the a;-space 
(ultra-)locality of our action is the same as for the unsmeared clover action - such a "normal" behavior is exactly 
what one would expect. To illustrate this point, a typical effective mass plot is shown in Fig.[6l Then, the masses are 
extracted from a correlated single channel cosh or sinh fit to the correlators. In order to estimate the systematic error 
due to excited states, we reduced the initial fit time by up to 2 timeslices and repeated the analysis with the new fit 
ranges. This difference then propagates into the systematic error in the continuum limit. 

For each coupling jS we then interpolate a^M^, aMp, aMjq and qMa linearly to a common current quark mass as 
determined by M^^/Mp. For illustration the interpolation at /? = 3.59 is shown in Fig. [71 The error on the current 
quark mass is of order 10~^ and therefore barely visible on this scale. Note that all data points are fully unquenched. 

We perform our scaling test on the baryon spectrum for three different values of M.^/Mp, all of which can be reached 
by interpolating our simulation data. In Tab.Uwe summarize the values of ampcAd oM^, aMp, oMn and aM^ after 
interpolation to M-j^jMp = 0.60, 0.64, 0.68. Also listed is LM^, which is roughly constant for fixed MT^/Mp. Moreover, 
even for the lightest data set we are deep in the M.^L > 4 regime. In case this criterion alone would not garantee the 
smallness of finite volume effects, the fact that our boxes have a fixed physical size ensures that such effects would be 
the same for all data at a given M^/Mp ratio, and the scaling test would still be meaningful. 

The masses are known to better than 2% and, due to correlations, this is also true for mass ratios. For the three 
lines of constant physics, Mn and Ma in units of Mtt are plotted in Fig. [8] as functions of the squared lattice spacing 
(see below), measured in units of the vector meson mass. We normalize the baryon masses by to clearly separate 
the lines of constant physics in the plot. The fits incorporate the error bars along both the vertical and horizontal 
axes. 

For both the spin-1/2 and spin-3/2 baryons, the continuum limit is approached smoothly with scaling violations 
of at most 1.2% at /? = 2.8. The extrapolations shown exclude this data point but consistent results are obtained by 
using all available data. 

While we expect that our choice of the clover coefficient is close to a non-perturbatively determined value, we cannot 
exclude effects that are linear in the lattice spacing in principle. The cutoff effects that we consider here are so small 
that we can not make a definitive statement, despite the fact that we have very precise data and cover more than a 
factor of seven in a^. Assuming the lattice artifacts to be linear in a results in an only marginally worse fit. 

An alternative way of proceeding is doing a combined chiral and continuum extrapolation with all datapoints 
at once. Applying this procedure one obtains basically consistent continuum limits and we assume the absolute 
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FIG. 7: Linear fits of tlie spectrum in terms of the PCAC quarlt mass. Data shown are from /3 = 3.59, L/a = 16 simulations. 
The line indicates the central value of the interpolation and the shaded region is the corresponding la error band (based on 
the assumption that the linear ansatz is correct). 
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0.0839(8) 5.03 
0.0607(23) 4.95 
0.0545(13) 5.12 
0.0405(6) 5.03 
0.0270(6) 5.41 


0.6292(21) 0.9832(33) 1.4341(43) 1.6581(59) 
0.4950(47) 0.7735(73) 1.127(10) 1.3074(82) 
0.4268(23) 0.6669(35) 0.9711(62) 1.1282(71) 
0.3146(23) 0.4916(36) 0.7099(35) 0.8278(28) 
0.2256(18) 0.3524(28) 0.5081(29) 0.5933(29) 


0.68 


8 2.80 
10 3.23 
12 3.40 
16 3.59 
24 3.76 


0.1050(11) 5.60 
0.0796(21) 5.57 
0.0693(12) 5.76 
0.0506(7) 5.57 
0.0343(9) 6.11 


0.6993(22) 1.0284(32) 1.5286(52) 1.7401(65) 
0.5574(52) 0.8198(76) 1.212(11) 1.389(10) 
0.4798(30) 0.7055(44) 1.0354(47) 1.1903(52) 
0.3483(22) 0.5122(32) 0.7495(30) 0.8590(41) 
0.2546(25) 0.3744(37) 0.5434(38) 0.6242(39) 



TABLE I: Results of the interpolation of aM,r, aMp, uMn and aM^, obtained from simulations performed at different bare 
quark masses and gauge couplings, to the reference points M^/Mp = 0.60,0.64,0.68. 

differences as our systematic errors. 

For illustrative purposes, we set the scale by linearly interpolating Mp and to the point where 

MjMp = ^2{Mf'^y - (AfP^y'')2/MP^'"' ^ 0.67 (7) 

and identify Mp with the mass of the physical 0. In this convention we cover lattice spacings from about 0.19 fm 
down to 0.07 fm (see Fig. IH]). In this range we find only small scaling violations in the spectrum and those disappear 
smoothly toward the continuum. The behavior is consistent with that of an 0(a)-improved theory. 
The scaling of other observables, especially matrix elements, will be investigated in the future. 

VII. SUMMARY 

We have described an efficient algorithm to perform full lattice QCD calculations with stout-link, improved clover 
fermions and demonstrated its potential with a scaling study of light baryon masses in Nf = 3 QCD. We have tested 
the algorithm and found it to be stable and reliable down to relatively coarse lattices with a ~ 0.16 fm. We have also 
monitored the stability of the MD integration and the lowest eigenvalue of the (even-odd preconditioned) fermion 
matrix and demonstrated that the latter is sufficiently far away from zero on all of our ensembles. Furthermore, 
we have shown that there is no sign of exceptional configurations even with substantially lighter pion masses in an 
AT^ = 2 -f 1 setting. 

Upon performing a "thermal cycle" at (3 — 3.3, with M^r ranging between ^ 240 MeV and ~ 440 MeV, we do not 
see any sign of a hysteresis. In other words, there is no indication of a nearby first order phase transition, even on 
fairly coarse lattices and for rather light quark masses. 

In a dedicated scaling test of light baryon masses, which included five lattice spacings with a total variation by 
almost a factor of three, we have demonstrated that scaling violations associated with the use of our stout-link clover 
action in full QCD are small for these quantities. Indeed, we have shown that discretization errors on light baryon 
masses do not exceed 2% for lattice spacings up to 0.19 fm. Moreover, all our data for a < 0.16 fm seem to be in the 
scaling window. This is in line with the findings of [33 | where a different approach to link smearing is taken. 

In conclusion, we find that the combination of a tree-level Symanzik improved gauge action and a six-step stout- 
smeared clover fermion action with csw = 1 is well suited for precision calculations of physical observables. We expect 
that the same will be true of other actions with comparable improvements. We look forward to presenting results for 
phenomenological quantities with this action in forthcoming papers. 
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FIG. 8: Mn and Ma, the mass of the spin-1/2 and spin-3/2 baryon, in terms of M^, versus the lattice spacing squared (in 
terms of Mp^). Each one of the three continuum extrapolations is based on the data at /3 = 3.76 — 3.23, but the curve is 
extended to /3 = 2.8 to allow for comparison. The continuum limits are Mn/M^ = 2.378(17)(43), 2.245(10)(51), 2.127(7)(34) 
and Ma/M^ = 2.827(23)(40), 2.626(17)(49), 2.446(16)(30) respectively. For all datapoints only statistical errors are shown. 
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FIG. 9: Scaling of the A and nucleon mass at Mj^/Mp — 0.67 in pliysical units using the scale setting procedure described 
around ©. In the continuum we obtain Mn = 1490(7)(27) MeV and Ma = 1720(10)(35) MeV. As in Fig.H only statistical 
errors are shown. 
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